pkgs <- c(
"here", "tidyverse", "RColorBrewer", "viridis", "brms", "patchwork", "tidylog",
"tidybayes", "lubridate", "tictoc", "modelr"
)
# minpack.lm needed if using nlsLM()
if (length(setdiff(pkgs, rownames(installed.packages()))) > 0) {
install.packages(setdiff(pkgs, rownames(installed.packages())), dependencies = T)
}
invisible(lapply(pkgs, library, character.only = T))
# Not on CRAN !
# devtools::install_github("seananderson/ggsidekick")
library(ggsidekick)
theme_set(theme_sleek())
# Set relative path
home <- here::here()Modelling temperature data as a covariate to length-at-age
Load packages
Temperature covariate
In this script, we fit smooth models to satellite data in order to predict average temperatures across years, lake stations during the growing season. We then combine these with older in-situ samples of lake temperature to create a full time series, which we then merge with the fish length data set.
Read satellite and in-situ data
# Load data
load(paste0(home, "/data/temp7_12.Rdata"))
head(temp_7to12)
#> # A tibble: 6 × 7
#> # Groups: monSiteCode, Lake_name, Date, year, yday [6]
#> monSiteCode Lake_name Date year yday ST satellite
#> <chr> <chr> <date> <dbl> <dbl> <dbl> <chr>
#> 1 LTL_Dr10 Druksiai_ne2 1984-04-12 1984 103 5.7 Landsat 5
#> 2 LTL_Dr10 Druksiai_ne2 1984-07-17 1984 199 16 Landsat 5
#> 3 LTL_Dr10 Druksiai_ne2 1984-10-12 1984 286 11 Landsat 5
#> 4 LTL_Dr10 Druksiai_ne2 1985-03-14 1985 73 0.7 Landsat 5
#> 5 LTL_Dr10 Druksiai_ne2 1985-05-17 1985 137 10.6 Landsat 5
#> 6 LTL_Dr10 Druksiai_ne2 1985-07-27 1985 208 18.8 Landsat 5
load(paste0(home, "/data/lst_final.RData"))
head(lst_final)
#> # A tibble: 6 × 7
#> # Groups: monSiteCode, Lake_name, Date, year, yday [6]
#> monSiteCode Lake_name Date year yday ST satellite
#> <chr> <chr> <date> <int> <int> <dbl> <chr>
#> 1 LTL_Dr0 Druksiai_warmmax 1984-04-12 1984 103 5.1 Landsat 5
#> 2 LTL_Dr0 Druksiai_warmmax 1984-05-21 1984 142 25.2 Landsat 5
#> 3 LTL_Dr0 Druksiai_warmmax 1985-08-12 1985 224 24.2 Landsat 5
#> 4 LTL_Dr0 Druksiai_warmmax 1985-08-28 1985 240 28.6 Landsat 5
#> 5 LTL_Dr0 Druksiai_warmmax 1986-03-17 1986 76 11.7 Landsat 5
#> 6 LTL_Dr0 Druksiai_warmmax 1986-05-04 1986 124 20 Landsat 5
d <- bind_rows(lst_final, temp_7to12) |>
ungroup() |>
mutate(lake_station = as.factor(Lake_name)) |>
drop_na(year) |>
rename(lst = ST)
# Clean covariates
d <- d |>
mutate(
yday = as.numeric(yday(Date)),
year = as.integer(year)
)
# Remove south-east bay (Druksiai_sebay)
unique(d$lake_station)
#> [1] Druksiai_warmmax Druksiai_warm Druksiai_ne Druksiai_nw
#> [5] Druksiai_west Druksiai_sebay Druksiai_ne2 Druksiai_west2
#> [9] Druksiai_sebay2 Druksiai_centre Druksiai_intake Druksiai_nw2
#> 12 Levels: Druksiai_centre Druksiai_intake Druksiai_ne ... Druksiai_west2
d |> filter(lake_station == "Druksiai_sebay")
#> # A tibble: 395 × 8
#> monSiteCode Lake_name Date year yday lst satellite lake_station
#> <chr> <chr> <date> <int> <dbl> <dbl> <chr> <fct>
#> 1 LTL_Dr5 Druksiai_seb… 1985-05-17 1985 137 12.9 Landsat 5 Druksiai_se…
#> 2 LTL_Dr5 Druksiai_seb… 1985-08-28 1985 240 19.9 Landsat 5 Druksiai_se…
#> 3 LTL_Dr5 Druksiai_seb… 1985-09-13 1985 256 13.5 Landsat 5 Druksiai_se…
#> 4 LTL_Dr5 Druksiai_seb… 1986-05-04 1986 124 13.3 Landsat 5 Druksiai_se…
#> 5 LTL_Dr5 Druksiai_seb… 1986-07-23 1986 204 22 Landsat 5 Druksiai_se…
#> 6 LTL_Dr5 Druksiai_seb… 1986-07-30 1986 211 24.6 Landsat 5 Druksiai_se…
#> 7 LTL_Dr5 Druksiai_seb… 1986-08-08 1986 220 22.8 Landsat 5 Druksiai_se…
#> 8 LTL_Dr5 Druksiai_seb… 1986-11-12 1986 316 3.6 Landsat 5 Druksiai_se…
#> 9 LTL_Dr5 Druksiai_seb… 1987-06-24 1987 175 22 Landsat 5 Druksiai_se…
#> 10 LTL_Dr5 Druksiai_seb… 1987-08-18 1987 230 15.4 Landsat 5 Druksiai_se…
#> # ℹ 385 more rows
d |> filter(lake_station == "Druksiai_sebay2")
#> # A tibble: 428 × 8
#> monSiteCode Lake_name Date year yday lst satellite lake_station
#> <chr> <chr> <date> <int> <dbl> <dbl> <chr> <fct>
#> 1 LTL_Dr12 Druksiai_seb… 1984-04-12 1984 103 5.1 Landsat 5 Druksiai_se…
#> 2 LTL_Dr12 Druksiai_seb… 1984-07-17 1984 199 16.7 Landsat 5 Druksiai_se…
#> 3 LTL_Dr12 Druksiai_seb… 1985-03-14 1985 73 3.7 Landsat 5 Druksiai_se…
#> 4 LTL_Dr12 Druksiai_seb… 1985-05-01 1985 121 6.4 Landsat 5 Druksiai_se…
#> 5 LTL_Dr12 Druksiai_seb… 1985-05-17 1985 137 14.6 Landsat 5 Druksiai_se…
#> 6 LTL_Dr12 Druksiai_seb… 1985-07-27 1985 208 19.4 Landsat 5 Druksiai_se…
#> 7 LTL_Dr12 Druksiai_seb… 1985-08-28 1985 240 22.3 Landsat 5 Druksiai_se…
#> 8 LTL_Dr12 Druksiai_seb… 1985-09-06 1985 249 21.9 Landsat 5 Druksiai_se…
#> 9 LTL_Dr12 Druksiai_seb… 1985-09-13 1985 256 18.5 Landsat 5 Druksiai_se…
#> 10 LTL_Dr12 Druksiai_seb… 1986-03-17 1986 76 1.4 Landsat 5 Druksiai_se…
#> # ℹ 418 more rows
d <- d |>
filter(!lake_station == "Druksiai_sebay") |>
droplevels()
# Insitu data
d_insit <- read.csv(paste0(home, "/data/Druksiai_temperature_raw.csv")) |>
filter(place == "station1") |>
mutate(
Date = paste(Year, Month, Day, sep = "-"),
Date = as.Date(Date),
yday = as.numeric(yday(Date)),
year = as.integer(Year)
) |>
rename("lst" = "Temp")Plot data
# Plot satellite data
ggplot(d, aes(yday, lst, color = lake_station)) +
geom_point() +
geom_smooth()
ggplot(d, aes(yday, lst, color = factor(year))) +
geom_point(alpha = 0.4) +
facet_wrap(~lake_station)
# Plot in-situ data
ggplot(d_insit, aes(yday, lst)) +
geom_point() +
geom_line() +
facet_wrap(~year)
ggplot(d_insit, aes(yday, lst, color = factor(year))) +
geom_point() +
geom_line()
ggplot(d, aes(yday)) +
geom_histogram() +
facet_wrap(~year)
# Combine
dd <- bind_rows(
d,
d_insit |>
dplyr::select(lst, yday, year, place) |>
rename(lake_station = place)
) |>
mutate(year_f = as.factor(year)) |>
dplyr::select(-monSiteCode, -Lake_name, -Date, -satellite)
mean_yday <- mean(dd$yday)
sd_yday <- sd(dd$yday)
unique(is.na(dd))
#> year yday lst lake_station year_f
#> [1,] FALSE FALSE FALSE FALSE FALSEFit satelite + in-situ model
In contrast to earlier versions, I pool in-situ and satelite. There isnt any overlap in stations, and stations differ in temperature, so it doesnt make sense to include “type” as a factor. I just treat it as a random site instead.
# s(yday): Seasonal pattern (shared)
# year_f: Population-average year effects
# (1 | lake_station): Persistent baseline difference for each station
# (1 | year_f:lake_station): Year-specific deviations by station
# m_prior <- brm(lst ~ s(yday, bs = "cc") +
# year_f +
# (1|lake_station) +
# (1|year_f:lake_station),
# data = dd,
# family = student(),
# iter = 1,
# chains = 1,
# knots = list(yday = c(0.5, 366.5)))
#
# get_prior(m_prior)
prior <- c(
prior(normal(0, 10), class = "b"),
prior(student_t(3, 18.2, 7.3), class = "Intercept"),
prior(student_t(3, 0, 7.3), class = "sds"),
prior(student_t(3, 0, 7.3), class = "sigma"),
prior(gamma(2, 0.1), class = "nu"),
prior(student_t(3, 0, 7.3), class = "sd", group = "lake_station"),
prior(student_t(3, 0, 7.3), class = "sd", group = "year_f:lake_station")
)
tic()
m <- brm(
lst ~ s(yday, bs = "cc") + year_f +
(1|lake_station) +
(1|year_f:lake_station),
data = dd,
prior = prior,
sample_prior = "yes",
family = student(),
iter = 5000,
chains = 4,
cores = 4,
knots = list(yday = c(0.5, 366.5)),
control = list(adapt_delta = 0.99, max_treedepth = 15)
)
#> Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
#> using C compiler: ‘Apple clang version 17.0.0 (clang-1700.6.3.2)’
#> using SDK: ‘MacOSX26.2.sdk’
#> clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/opt/R/arm64/include -fPIC -falign-functions=64 -Wall -g -O2 -c foo.c -o foo.o
#> In file included from <built-in>:1:
#> In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
#> In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
#> In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
#> /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
#> 679 | #include <cmath>
#> | ^~~~~~~
#> 1 error generated.
#> make: *** [foo.o] Error 1
toc()
#> 1509.001 sec elapsed
# ~20 minutesCheck model fit
pp_check(m) +
scale_color_manual(
values = viridis(n = 5, option = "plasma")[c(2, 4)],
labels = c("obs", "yrep")
) +
theme(
legend.position.inside = c(0.05, 0.95),
legend.title = element_blank()
) +
labs(x = "Value", y = "Density") +
theme_sleek() +
guides(color = guide_legend(position = "inside")) +
theme(
legend.position.inside = c(0.1, 0.93),
legend.title = element_blank()
) +
NULL
ggsave(paste0(home, "/figures/supp/pp_check_temp.png"), width = 11, height = 11, units = "cm")Plot prior vs posterior
# Get all posteriors (remove helper columns and prior columns)
draws <- as_draws_df(m)
posteriors <- draws |>
dplyr::select(
-starts_with("prior_"),
-lprior, -lp__,
-.chain, -.iteration, -.draw,
-starts_with("r_f_month"),
-Intercept
) |>
pivot_longer(everything(),
names_to = "parameter",
values_to = "value"
) |>
mutate(
parameter = case_when(
parameter == "sd_lake_station__Intercept" ~ "sd_lake_station",
parameter == "sd_year_f:lake_station__Intercept" ~ "sd_year_f:lake_station",
TRUE ~ parameter
),
type = "posterior"
)
#> Warning: Dropping 'draws_df' class as required metadata was removed.
#> pivot_longer: reorganized (b_Intercept, b_year_f1979, b_year_f1980, b_year_f1981, b_year_f1982, …) into (parameter, value) [was 10000x518, now 5180000x2]
#> mutate: changed 20,000 values (<1%) of 'parameter' (0 new NA)
#> new variable 'type' (character) with one unique value and 0% NA
# Get all priors
priors <- draws |>
dplyr::select(starts_with("prior_"), -lprior) |>
pivot_longer(everything(),
names_to = "parameter",
values_to = "value"
) |>
mutate(
parameter = str_remove(parameter, "^prior_"),
parameter = case_when(
parameter == "Intercept" ~ "b_Intercept",
TRUE ~ parameter
),
type = "prior"
)
#> Warning: Dropping 'draws_df' class as required metadata was removed.
#> pivot_longer: reorganized (prior_Intercept, prior_b, prior_sds_syday, prior_sigma, prior_nu, …) into (parameter, value) [was 10000x7, now 70000x2]
#> mutate: changed 70,000 values (100%) of 'parameter' (0 new NA)
#> new variable 'type' (character) with one unique value and 0% NA
prior_post <- bind_rows(posteriors, priors)
# Group parameters to plot together
prior_post_clean <- prior_post |>
mutate(
par_group = case_when(
str_starts(parameter, "b_Intercept") ~ "intercept",
str_starts(parameter, "b_y") ~ "b",
str_starts(parameter, "b") ~ "b",
# str_starts(parameter, "s_") ~ "s_sds",
# str_starts(parameter, "sds_") ~ "s_sds",
str_starts(parameter, "sigma") ~ "sigma_sd_nu",
str_starts(parameter, "sd_") ~ "sigma_sd_nu",
str_starts(parameter, "nu") ~ "sigma_sd_nu",
TRUE ~ NA
)
) |>
drop_na(par_group)
#> mutate: new variable 'par_group' (character) with 4 unique values and 90% NA
#> drop_na: removed 4,700,000 rows (90%), 550,000 rows remaining
for (i in unique(prior_post_clean$par_group)) {
d_sub <- prior_post_clean |>
filter(par_group == i)
limits <- d_sub |>
filter(type == "posterior") |>
group_by(parameter) |>
summarise(
lower = quantile(value, 0.05),
upper = quantile(value, 0.975)
)
ncol <- if (i == "sigma_sd_nu") {
2
} else if (i == "b") {
10
} else {
4
}
if (i == "b") {
d_sub |>
left_join(limits, by = "parameter") |>
# filter(value > 0.1 * lower, value < upper * 10) |>
ggplot(aes(value)) +
geom_density(
data = \(d) d |> filter(type == "posterior"),
aes(fill = parameter, color = parameter),
alpha = 0.4,
color = NA
) +
geom_density(
data = \(d) d |> filter(type == "prior"),
aes(group = parameter),
fill = "grey60",
alpha = 0.4,
color = "grey20"
) +
scale_fill_viridis(option = "plasma", discrete = TRUE) +
scale_color_viridis(option = "plasma", discrete = TRUE) +
labs(x = "Value", y = "Density") +
coord_cartesian(expand = 0) +
theme(
strip.text.x.top = element_text(size = 8),
legend.position = "bottom",
legend.title = element_blank()
)
} else {
(
d_sub |>
left_join(limits, by = "parameter") |>
filter(value > 0.1 * lower, value < upper * 10) |>
ggplot(aes(value)) +
geom_density(
data = \(d) d |> filter(type == "posterior"),
aes(fill = parameter, color = parameter),
alpha = 0.4,
color = NA
) +
geom_density(
data = \(d) d |> filter(type == "prior"),
aes(group = parameter),
fill = "grey60",
alpha = 0.4,
color = "grey20"
) +
scale_fill_viridis(option = "plasma", discrete = TRUE) +
scale_color_viridis(option = "plasma", discrete = TRUE) +
labs(x = "Value", y = "Density") +
coord_cartesian(expand = 0) +
facet_wrap(~parameter, ncol = ncol, scales = "free") +
theme(
strip.text.x.top = element_text(size = 8),
legend.position = "bottom",
legend.title = element_blank()
)
)
}
ggsave(paste0(home, "/figures/supp/prior_post/prior_post_temp_", i, ".png"), width = 18, height = 22, units = "cm")
}
#> filter: removed 530,000 rows (96%), 20,000 rows remaining
#> filter: removed 10,000 rows (50%), 10,000 rows remaining
#> group_by: one grouping variable (parameter)
#> summarise: now one row and 3 columns, ungrouped
#> left_join: added 2 columns (lower, upper)
#> > rows only in x 0
#> > rows only in y ( 0)
#> > matched rows 20,000
#> > ========
#> > rows total 20,000
#> filter: removed 510 rows (3%), 19,490 rows remaining
#> filter: removed 9,490 rows (49%), 10,000 rows remaining
#> filter: removed 10,000 rows (51%), 9,490 rows remaining
#> filter: removed 100,000 rows (18%), 450,000 rows remaining
#> filter: removed 10,000 rows (2%), 440,000 rows remaining
#> group_by: one grouping variable (parameter)
#> summarise: now 44 rows and 3 columns, ungrouped
#> left_join: added 2 columns (lower, upper)
#> > rows only in x 10,000
#> > rows only in y ( 0)
#> > matched rows 440,000
#> > =========
#> > rows total 450,000
#> filter: removed 10,000 rows (2%), 440,000 rows remaining
#> filter: removed 440,000 rows (98%), 10,000 rows remaining
#> filter: removed 470,000 rows (85%), 80,000 rows remaining
#> filter: removed 40,000 rows (50%), 40,000 rows remaining
#> group_by: one grouping variable (parameter)
#> summarise: now 4 rows and 3 columns, ungrouped
#> left_join: added 2 columns (lower, upper)
#> > rows only in x 0
#> > rows only in y ( 0)
#> > matched rows 80,000
#> > ========
#> > rows total 80,000
#> filter: removed 3,385 rows (4%), 76,615 rows remaining
#> filter: removed 36,615 rows (48%), 40,000 rows remaining
#> filter: removed 40,000 rows (52%), 36,615 rows remainingPredict from model
# Annual average growth season temperature by station
gs_min <- yday("2020-05-01")
gs_max <- yday("2020-10-15")
# This is a pretty big prediction grid, so therefore I need to do it by site, then summarise, else the data frame explodes in size
nd <- dd |>
distinct(lake_station, year_f) |>
expand_grid(yday = 1:365)
pred_list <- list()
for (i in unique(dd$lake_station)) {
nd_sub <- nd |> filter(lake_station == i)
pred_list[[i]] <- nd_sub |>
add_epred_draws(m) |>
group_by(lake_station, year_f, yday) |>
median_qi(.epred, .width = c(0.8, 0.95)) |>
ungroup()
}
pred_sum <- bind_rows(pred_list) |>
mutate(year = as.numeric(as.character(year_f)))
pred_sum |>
filter(.width == 0.95) |> # this is just because the tibble is duplicated for each width
ggplot(aes(yday, .epred, color = lake_station)) +
geom_line() +
facet_wrap(~year) +
geom_point(data = d, aes(yday, lst, color = lake_station)) +
theme(legend.position = "bottom") +
scale_color_viridis(name = "Lake station", discrete = TRUE)
# Heatmap
pred_sum |>
filter(.width == 0.95) |>
ggplot(aes(yday, y = year, fill = .epred, color = .epred)) +
geom_raster() +
facet_wrap(~lake_station, ncol = 3) +
scale_fill_viridis(option = "magma") +
scale_color_viridis(option = "magma") +
labs(x = "Day of the year", y = "Year", color = "Predicted LST (°C)", fill = "Predicted LST (°C)") +
coord_cartesian(expand = 0) +
guides(fill = guide_colorbar(
title.position = "top", title.hjust = 0.5
)) +
theme(
legend.key.width = unit(0.9, "cm"),
legend.position = "bottom",
)
ggsave(paste0(home, "/figures/supp/temp_station_heat.png"), width = 17, height = 22, units = "cm")
# Summarise temp for a given time period and plot
d_sum <- dd |>
filter(yday > gs_min & yday < gs_max) |>
group_by(year, lake_station) |>
summarise(mean_temp = mean(lst)) |>
ungroup()
sum_pred_g <- pred_sum |>
# filter(.width == 0.95) |>
filter(yday > gs_min & yday < gs_max) |>
group_by(year, lake_station, .width) |>
summarise(
lwr = mean(.lower),
est = mean(.epred),
upr = mean(.upper)
) |>
ungroup()
sum_pred_g |>
mutate(lake_station = str_remove(lake_station, "Druksiai_")) |>
ggplot(aes(year, est,
linetype = "Model (95% C.I.)",
color = lake_station, fill = lake_station
)) +
facet_wrap(~lake_station, ncol = 3) +
geom_ribbon(
data = \(d) d |> filter(.width == 0.95),
aes(ymin = lwr, ymax = upr),
color = NA, alpha = 0.3
) +
# geom_ribbon(data = \(d) d |> filter(.width == 0.8),
# aes(ymin = lwr, ymax = upr),
# color = NA, alpha = 0.1) +
geom_point(
data = d_sum |> mutate(lake_station = str_remove(lake_station, "Druksiai_")),
aes(year, mean_temp, shape = "Data"),
color = "grey30", inherit.aes = FALSE, size = 0.8
) +
geom_line() +
scale_fill_viridis(option = "plasma", discrete = TRUE) +
scale_color_viridis(option = "plasma", discrete = TRUE) +
guides(
shape = guide_legend(
override.aes = list(size = 2)
),
fill = "none", color = "none"
) +
scale_y_continuous(breaks = seq(5, 40, by = 2)) +
labs(
x = "Year",
y = expression("Mean lake surface temperature (°C) (May 1"^"st" * "–Oct 15"^"th" * ")")
) +
theme(
legend.position = "bottom",
legend.spacing.y = unit(-0.2, "cm"),
legend.title = element_blank()
)
ggsave(paste0(home, "/figures/supp/temp_pred.png"), width = 17, height = 19, units = "cm")
# Make a polar plot?
pred_sum |>
filter(.width == 0.95) |>
summarise(.epred = mean(.epred), .by = c(lake_station, yday)) |>
summary(.epred)
#> lake_station yday .epred
#> Length:4380 Min. : 1 Min. :-1.885
#> Class :character 1st Qu.: 92 1st Qu.: 2.263
#> Mode :character Median :183 Median : 9.554
#> Mean :183 Mean :10.658
#> 3rd Qu.:274 3rd Qu.:19.445
#> Max. :365 Max. :27.644
pred_sum |>
filter(.width == 0.95) |>
summarise(.epred = mean(.epred), .by = c(lake_station, yday)) |>
ggplot(aes(yday, .epred, color = lake_station)) +
geom_line(linewidth = 1) +
coord_polar(start = 0, direction = -1) +
scale_x_continuous(
breaks = c(1, 32, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335),
labels = c("Jan", "Feb", "Mar", "Apr", "May", "Jun",
"Jul", "Aug", "Sep", "Oct", "Nov", "Dec"),
limits = c(0, 365)
) +
scale_y_continuous(
breaks = seq(0, 30, by = 5)
) +
annotate("text", x = 1, y = 5, label = "5°C", size = 3.5, color = "grey40") +
annotate("text", x = 1, y = 10, label = "10°C", size = 3.5, color = "grey40") +
annotate("text", x = 1, y = 15, label = "15°C", size = 3.5, color = "grey40") +
annotate("text", x = 1, y = 20, label = "20°C", size = 3.5, color = "grey40") +
annotate("text", x = 1, y = 25, label = "25°C", size = 3.5, color = "grey40") +
scale_color_viridis_d(option = "magma", end = 0.9) +
theme_minimal() +
theme(
axis.text.y = element_blank(),
axis.title.y = element_blank(),
axis.title.x = element_blank()
) +
labs(
color = "Lake Station"
)Now combine the series into a single that that we use a co-variate.
annual_trends <- dd |>
distinct(year) |>
mutate(year_f = as.factor(year)) |>
expand_grid(
yday = gs_min:gs_max
) |>
add_epred_draws(
m,
re_formula = NA # average site
) |>
ungroup() |>
# Summarise across days-of-the day by averaging by earch draw and year only
summarise(
mean_lst = mean(.epred),
.by = c(.draw, year_f)
) |>
# Summarise across draws
group_by(year_f) |>
median_qi(mean_lst, .width = c(0.8, 0.95)) |>
ungroup() |>
mutate(year = as.numeric(as.character(year_f)))
ggplot(annual_trends, aes(year, mean_lst)) +
geom_vline(xintercept = 1981, linetype = 2, alpha = 0.5) +
geom_vline(xintercept = 2009, linetype = 2, alpha = 0.5) +
geom_ribbon(
data = \(d) d |> filter(.width == 0.8),
aes(ymin = .lower, ymax = .upper),
color = NA, alpha = 0.3, fill = "grey50"
) +
geom_ribbon(
data = \(d) d |> filter(.width == 0.95),
aes(ymin = .lower, ymax = .upper),
color = NA, alpha = 0.3, fill = "grey70"
) +
geom_line() +
scale_y_continuous(breaks = seq(5, 40, by = 1)) +
labs(x = "Year", y = expression("Mean lake surface temperature (°C) (May 1"^"st" * "–Oct 15"^"th" * ")")) +
coord_cartesian(expand = 0) +
annotate("text", x = 1979.5, y = 23, label = "Pre\nwarming\n(a)", size = 2.5) +
annotate("text", x = 1991, y = 23, label = "During operation of power plant (b)", size = 2.5) +
annotate("text", x = 2016, y = 23, label = "Post operation of power plant (c)", size = 2.5) +
annotate(
"rect",
xmin = 1988, xmax = 2004,
ymin = -Inf, ymax = Inf,
alpha = 0.08,
fill = "grey50"
) +
annotate(
"text",
x = 1996, y = 13.6,
label = "2 reactors",
size = 2.5,
fontface = "italic"
) +
annotate(
"text",
x = 1984.5, y = 13.6,
label = "1 reactor",
size = 2.5,
fontface = "italic"
) +
annotate(
"text",
x = 2006.5, y = 13.6,
label = "1 reactor",
size = 2.5,
fontface = "italic"
) # annotate(
# "segment",
# x = 1988, xend = 1988,
# y = 14.1, yend = 13.3,
# arrow = arrow(length = unit(0.15, "cm")),
# linewidth = 0.4
# ) +
# annotate(
# "text",
# x = 1988, y = 14.3,
# label = "Second reactor initiated",
# size = 2.5,
# hjust = 0.5
# ) +
# annotate(
# "segment",
# x = 2004, xend = 2004,
# y = 14.1, yend = 13.3,
# arrow = arrow(length = unit(0.15, "cm")),
# linewidth = 0.4
# ) +
# annotate(
# "text",
# x = 2004, y = 14.3,
# label = "Second reactor closed",
# size = 2.5,
# hjust = 0.5
# )
ggsave(paste0(home, "/figures/temp_series.png"), width = 17, height = 11, units = "cm")write_csv(annual_trends, paste0(home, "/output/merged_temp_covar_data.csv"))